clear all;
clc;
tstart=1; %start time
tmax=50; %max time of simulation
time=tstart:tmax; %time vector

n_iterations=5;
%nparticles = 1000;
nparticles = 10;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(npos,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(npos,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(npos,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(npos,npos);%a's for spin down at each N at current time

% T=tmax/2:100:tmax;
T=time;
pos_n_counter=0;
pos_m_counter=0;
a=0;
epsilon = 0.025;
p_n_tmax = zeros(npos,npos);
p_n_u_tmax = zeros(npos,npos);
p_n_d_tmax = zeros(npos,npos);

p_n = zeros(npos,npos);
p_n_u= zeros(npos,npos);
p_n_d= zeros(npos,npos);

for iter=1:n_iterations
    display(iter);
    a = 1/(sqrt(2));
%     a = rand;
    a_n_u_prev(nmid_rel,nmid_rel)=a;
%     a_n_d_prev(nmid_rel,nmid_rel)=(sqrt(1-a^2))*i;
    a_n_d_prev(nmid_rel,nmid_rel)=a*i;
    
    
    %for all times
    for t=time
        display(t);
        pos_m_counter = 0;
        pos_n_counter = 0;
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
%         if(t ~= 1)
% %              a=1/sqrt(2);
%              a=(1/(sqrt(2))) + epsilon*rand;
% %             a = epsilon*rand;
%         end;

        %for all possible positions for this particular time, fill up the a
        %arrays for the current time
            
        for pos_n_counter=nstart_rel:nend_rel
%             display(pos_n_counter);
            for pos_m_counter=nstart_rel:nend_rel
%                   display(pos_m_counter);
                if(rand < 0.5)
                    a_n_u_cur(pos_n_counter,pos_m_counter)= 0.5*(a_n_u_prev(pos_n_counter-1,pos_m_counter) *(a + sqrt(1-a^2)) + a_n_d_prev(pos_n_counter-1,pos_m_counter) * (-a + sqrt(1-a^2)) + a_n_d_prev(pos_n_counter+1,pos_m_counter) * (a-sqrt(1-a^2)) + a_n_d_prev(pos_n_counter + 1,pos_m_counter) * (a+sqrt(1-a^2)) );
                    a_n_d_cur(pos_n_counter,pos_m_counter) = 0.5*(a_n_u_prev(pos_n_counter-1,pos_m_counter) * (a+sqrt(1-a^2)) + a_n_d_prev(pos_n_counter-1,pos_m_counter)*(-a + sqrt(1-a^2)) + a_n_u_prev(pos_n_counter+1,pos_m_counter)*(-a + sqrt(1-a^2)) - a_n_d_prev(pos_n_counter + 1,pos_m_counter) * (a+sqrt(1-a^2)) );
                    
%                     a_n_u_cur(pos_n_counter,pos_m_counter)=  ( sqrt(1-a*a) * a_n_d_prev(pos_n_counter-1,pos_m_counter) )+ ( a * a_n_u_prev(pos_n_counter-1,pos_m_counter) );
%                     a_n_d_cur(pos_n_counter,pos_m_counter)=( sqrt(1-a*a ) * a_n_u_prev(pos_n_counter+1,pos_m_counter)) - ( a * a_n_d_prev(pos_n_counter+1,pos_m_counter) );
                    
                else
                    a_n_u_cur(pos_n_counter,pos_m_counter)= 0.5*(a_n_u_prev(pos_n_counter,pos_m_counter-1) *(a + sqrt(1-a^2)) + a_n_d_prev(pos_n_counter,pos_m_counter-1) * (-a + sqrt(1-a^2)) + a_n_d_prev(pos_n_counter,pos_m_counter+1) * (a-sqrt(1-a^2)) + a_n_d_prev(pos_n_counter,pos_m_counter+1) * (a+sqrt(1-a^2)) );                    
                    a_n_d_cur(pos_n_counter,pos_m_counter) = i*0.5*(a_n_u_prev(pos_n_counter,pos_m_counter-1) * (a+sqrt(1-a^2)) + a_n_d_prev(pos_n_counter,pos_m_counter-1)*(-a + sqrt(1-a^2)) + a_n_u_prev(pos_n_counter,pos_m_counter+1)*(-a + sqrt(1-a^2)) - a_n_d_prev(pos_n_counter,pos_m_counter+1) * (a+sqrt(1-a^2)) );                    
                end;

                %record the probability densities of up,down and all particles
%                 p_n_u_tmax(pos_n_counter,pos_m_counter) = (abs(a_n_u_cur(pos_n_counter,pos_m_counter)))^2;
%                 p_n_d_tmax(pos_n_counter,pos_m_counter) = (abs(a_n_d_cur(pos_n_counter,pos_m_counter)))^2; 
%                 p_n_tmax(pos_n_counter,pos_m_counter) = (abs(a_n_u_cur(pos_n_counter,pos_m_counter)))^2 + (abs(a_n_d_cur(pos_n_counter,pos_m_counter)))^2; 
            end;
        end;
        
        %record the probability densities of up,down and all particles
        p_n_u_tmax = abs(a_n_u_cur).^2;
        p_n_d_tmax = abs(a_n_d_cur).^2; 
        p_n_tmax =  p_n_u_tmax + p_n_d_tmax;     
        
        p_n = p_n + p_n_tmax;
        p_n_u = p_n_u + p_n_u_tmax;
        p_n_d = p_n_d + p_n_d_tmax;
        
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
        if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(npos,npos);
            a_n_d_cur = zeros(npos,npos);
        end;
        
        if(iter == n_iterations-1) 
            p_n_tmax = zeros(npos,npos);
            p_n_u_tmax = zeros(npos,npos);
            p_n_d_tmax = zeros(npos,npos);
        end;
        
    end;
   
end;

p_n = p_n/(n_iterations*tmax);
p_n_u = p_n_u/(n_iterations*tmax);
p_n_d = p_n_d/(n_iterations*tmax);
% 

figure(10)
surf(pos-nmid_rel,pos-nmid_rel,p_n);
title('plot of P_n(t) for 2D Quantum hadamard random walk');
grid on;

figure(11)
surf(pos-nmid_rel,pos-nmid_rel,p_n_u);
title('plot of P_n(t) for 2D Quantum hadamard random walk - up-spin');
grid on;